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For  a realistic  and  practical  aerodynamic  optimization  the  most  appropriate  combination  of  the  three  sets  of  tools 
taking  part  in  the  process  should  be  carefully  studied.  That  is,  the  optimization  should  allow  an  easy  implementation 
of  constraints,  and  should  be  multipoint  without  the  need  to  prescribe  pressure  distributions  in  the  objective  function; 
the  design  space  should  be  broad  enough;  and  the  analysis  tool  should  be  fast  and  robust.  Taking  into  account  these 
criteria,  a code  for  multipoint  design  and  optimization  of  wings  in  subsonic  and  transonic  regime  has  been  developed 
and  will  be  described  in  this  paper.  The  objective  can  be  any  combination  of  the  global  aerodynamic  coefficients, 
and  geometrical  and  physical  constraints  can  be  applied.  Results  for  subsonic  and  transonic  cases  will  be  presented. 
Flexibility  in  the  use  of  the  design  variables  allows  many  different  tests  to  be  performed  before  the  best  solution  is 
achieved.  Lastly,  the  computational  cost  is  reduced  by  the  use  of  a low  level  code  for  computing  the  aerodynamic 
coefficients. 


1 Introduction 

The  quality  of  the  results  in  the  use  of  numerical  opti- 
mization for  the  design  of  aerodynamic  configurations 
depends  on  the  appropriate  choice  of  a lot  of  parameters. 
These  are  usually  interdependent  and  do  not  always 
affect  the  solution  in  the  same  direction,  especially  when 
the  gradient  computation  is  required.  For  example,  the 
use  of  more  complex  models  for  computing  the  flow 
does  not  necessarily  improve  the  results  because  of  the 
numerical  noise;  this  also  applies  to  the  type  and  number 
of  design  variables;  as  well  as  the  way  in  which  the 
global  aerodynamic  coefficients  are  determined,  whether 
by  local  or  far  field  methods,  among  many  other  factors. 

The  first  decision  to  be  taken  concerns  the  type  and 
number  of  design  variables,  which  is  a subject  that  has 
not  been  sufficiently  covered  in  our  opinion.  The  design 
space  should  be  broad  enough  to  avoid  the  dependency 
of  the  starting  geometry  in  the  final  solution  as  much 
as  possible.  As  a first  step  we  have  found  it  is  very 
important  to  replace  the  set  of  points  usually  defining 
the  initial  geometry  by  an  analytical  definition.  For  that 
purpose,  an  automatic  adjustment  by  Bezier  polynomia 
has  been  used.  After  the  validation  with  a data  base  of 
many  different  airfoils  it  was  concluded  that  at  least  25 
control  points  are  required  to  have  an  an  error  of  defini- 
tion less  than  0.01%,  which  seems  a reasonable  criteria  if 
the  recommended  initial  changes  in  the  design  variables 
are  0.001.  However,  the  shape  functions  originally 
included  in  the  code  (Wagner,  Hicks- Vanderplaats,  Leg- 
endre, etc  . . .)  are  still  available  for  comparison  purposes. 

As  regards  the  CFD  codes  to  be  used  for  computing 
the  objective  function  robustness  is  their  most  essential 
characteristic  for  using  them  in  an  optimization  loop,  so 
that  their  sensitivity  to  the  numerical  data  can  be  kept 
as  low  as  possible.  For  this  reason,  and  not  only  for  the 
computing  time,  it  is  convenient  to  use  low  level  codes 
as  long  as  the  nature  of  the  flow  allows  us  to  do  so,  and 
particularly  if  the  objective  is  a function  of  the  global 
aerodynamic  coefficients  and,  almost  in  any  case,  when 
we  are  in  the  first  stages  of  the  optimization,  far  from 
the  final  solution. 


For  the  optimization  itself  a finite  differences  based 
commercial  code  [1]  has  been  used  with  the  options  of 
following  the  steepest  descent  or  conjugate  gradients  to 
where  we  have  added  the  option  of  using  a quasi-Newton 
method  [2]. 

These  three  sets  of  numerical  tools  have  been  combined 
in  a modular  way  to  develop  an  optimization  code  for 
airfoils  with  high  lift  devices  ( OPTPER ) and  another 
one  for  wings  ( OPT  WIN ) that  were  employed  for  the 
application  cases  to  be  described  in  this  report.  They 
allow  us  to  choose  among  different  CFD  codes  ranging 
from  velocity  potential  coupled  with  boundary  layer 
until  Navier- Stokes  (2D)  and  until  full  potential  coupled 
with  boundary  layer  (3D).  The  objective  function  can  be 
any  combination  of  the  global  aerodynamic  coefficients, 
and  the  equality  or  inequality  restrictions  to  be  imposed 
can  be  geometrical  (leading  edge  radius,  trailing  edge 
angle,  maximum  thickness,  camber,  etc. . .)  or  physical 
(pitching  moment  coefficient,  minimum  pressure,  etc. . .). 
Also  we  profited  from  the  expertise  collected  in  data 
bases  such  as  ESDU  for  checking  some  other  con- 
straints like  the  maximum  lift  coefficient  or  the  highest 
adverse  pressure  gradient  along  the  optimization  process. 

2 Optimization  procedure 

The  main  code  has  been  developed  in  a modular  form  in 
order  to  allow  easy  improvement  or  modification  of  any 
of  its  elements. 

In  general  a geometry  exists  that  requires  aerodynamical 
improvement  at  certain  flow  conditions,  but  it  is  also 
possible  to  start  the  computation  from  any  wing  section 
selected  from  a given  data  base  of  airfoils,  which  always 
gives  a kind  of  a guarantee  that  the  fin?il  solution  does 
not  depend  much  on  the  initial  wing.  Its  planform 
can  also  be  optimized  during  the  minimization  process, 
but  it  is  mainly  designed  using  structural  or  handling 
qualities  reasons. 

Next  the  objective  function  has  to  be  specified.  It  can  be 
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Figure  1:  Flow  chart  for  the  aerodynamic  optimiza- 
tion procedure 


new  geometry  is  recorded  at  the  end  of  every  iteration, 
the  user  may  stop  the  computation  at  any  time  without 
loosing  the  previous  results,  from  which  a new  design 
strategy  may  start.  In  order  to  verify  the  quality  of 
the  final  design  a drag  polar  is  computed  at  the  end, 
and  compared  with  that  of  the  original  geometry.  The 
procedure  described  above  is  illustrated  in  figure  1. 

3 Design  variables 

Concerning  the  design  variables  different  possibilities 
have  been  included.  First,  a set  of  shape  functions  to 
be  added  to  the  original  wing  section  was  used:  Wag- 
ner, Hicks- Vanderplaats,  Legendre,  cubic  patches,  etc  . . . 
After  several  trials  the  following  functions  were  chosen 
[3]: 

• Thickness 


ft(x)  = 4 • (l  - (1) 

being  p(x  m,  ) 

log 

P(Xmi)  = log(*m<)  0 " ~ 1 (2) 

and  Xmi  is  the  percentage  of  chord  where  the  max- 
imum of  ft(x)  is  located. 

• Camber 


any  combination  of  the  global  aerodynamic  coefficients 
at  one  or  several  points  of  the  drag  polar.  Apart 
from  this,  a pressure  distribution  can  be  prescribed  as 
objective,  but  for  this  purpose  the  inverse  methods  are 
generally  more  efficient. 

Caryying  out  an  exclusively  aerodynamic  optimization 
other  disciplines  are  not  taken  into  account  (structures, 
etc...).  Generally  it  is  also  necessary  to  impose  geo- 
metrical and/or  physical  constraints.  At  present  it  is 
possible  to  prescribe  limits  in  the  wing  section  area, 
leading  edge  radius,  angles  of  thickness  and  camber  line 
at  the  trailing  edge,  thickness  at  a given  positions  along 
the  chord,  minimum  lift,  maximum  drag,  absolute  value 
of  pitching  moment  and  maximum  velocity.  This  is  also 
important  in  order  to  avoid  unfeasible  geometries,  which 
is  always  a risk  when  the  design  space  is  sufficiently  large. 

For  computing  the  objective  function  several  analysis 
codes  can  be  selected.  In  this  report  a 3D  full  potential 
method  coupled  with  an  integral  boundary  layer  method 
has  been  used.  It  will  be  described  in  the  following 
chapters  along  with  the  design  variables. 

As  a minimization  tool  the  original  CONMIN  code  cou- 
pled with  a BFGS  routine  to  accelerate  the  convergence 
has  been  implemented  (CONMIN2). 

As  regard  as  the  data  input,  it  was  distributed  in  the 
following  way:  flow  conditions  and  objective  values, 

initial  geometry,  numerical  parameters  of  the  analysis, 
numerical  parameters  of  the  optimization,  and  number 
and  type  of  design  variables. 

Before  the  minimization  process  starts  the  user  may 
improve  the  initial  geometry  by  distributing,  smoothing 
and  interpolating  the  original  points  if  there  are  nu- 
merical oscillations  in  the  curvature,  which  is  internally 
evaluated. 


/«(*)  = 4 • (1  - z)p(*mi)  (l  - (1  - (3) 


being  p(®mi) 


P{xmi  ) 


lpg  (I) 

log  (1  - smi) 


0 < xmi  < i 


(4) 


whereas  for  | < xmi  < 1 the  expressions  1 and  2 
are  used. 


They  provide  a more  general  design  space  than  the  other 
functions,  and  the  number  of  design  variables  can  be 
easily  increased  by  an  approppriate  election  of  £mi.  An 
infinite  slope  at  the  leading  edge  for  the  camber  line  has 
been  avoided  (fig.  2). 

In  this  case  we  have  considered  nine  design  variables  in 
order  to  control  the  radius  of  the  leading  edge,  position 
and  value  of  the  maximum  thickness  and  camber,  angle 
of  thickness  and  camber  at  the  trailing  edge,  slope  of 
the  camber  line  at  the  leading  edge  and  with  at  least 
one  more  degree  of  freedom. 

However,  when  described  as  usual  by  a set  of  coordinates 
the  numerical  uncertainties  in  the  definition  of  the  origi- 
nal airfoil  can  have  a negative  effect  in  the  optimization. 
It  seems  convenient  to  have  the  airfoil  itself  defined  from 
the  begining  by  analytical  functions,  and  not  only  by  the 
deformation  shape  functions.  To  this  end  a method  for 
fitting  the  airfoil  to  a set  of  Bezier  polynomia  has  been 
developed  and  implemented  into  the  code  [4].  The  con- 
trol points  are  then  used  as  design  variables. 

One  function  is  used  for  each  upper  and  lower  side  of  the 
airfoil,  where  at  the  leading  edge  the  continuity  of  first 
and  second  derivatives  is  imposed: 


Lastly,  the  following  information  can  be  obtained  from 
the  output:  final  geometry,  evolution  of  design  variables, 
aerodynamic  coefficients  and  objective  function.  As  the 


NB 

v».i(«)  = ( Nf  ) **(i  - °)NB- 


(5) 
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Figure  2:  Thickness  and  camber  shape  functions 


being  Y the  fitted  coordinates,  NB  the  number  of  Bezier 
control  points,  s ranges  between  0, . . . , 1.  The  subscripts 
u and  l denote  upper  and  lower  side  respectively,  whilst 
the  superscript  B means  the  Bezier  control  point  coordi- 
nates. 

By  solving  two  least  square  problems  the  control  points 
Y*i  are  computed.  This  can  be  done  by  numerical  op- 
timization or  by  directly  solving  the  system  of  algebraic 
equations: 

[Aij]{YB}  = {BYi}  = — l (6) 

where, 

= £ [( N,s ) ( Nf ) •£”<*  - ..r-'H 

fc= 1 


NP 

Bn  = ^ [(  ) 4(1  - sk)NB~'  ■ (yo{sk)~ 

fc=l 

j/o(Si)(l-S*)WB-yo(^p)^s)]  (8) 

NP  is  the  number  of  points  defining  the  wing  section 
and  yo  are  the  coordinates  of  the  original  wing  sections. 

We  have  used  this  last  approach  not  only  because  of 
the  much  lower  computing  time  required  but  also  and 
mainly  because  not  so  many  numerical  parameters  are 
required  as  in  the  first  method. 

Once  the  problem  is  posed  and  the  solution  method 
developed  the  main  question  is  about  the  most  appro- 
priate number  of  control  points  for  defining  the  airfoil 
and  its  deformations  in  order  to  have  a large  enough 
design  space.  A set  of  very  different  airfoils  has  been 
tested  (conventional,  laminar,  supercritical,  high-lift, 
etc. . .)  and  the  outcome  is  that  if  the  average  error 
needs  to  be  less  than  10~4,  around  25  control  points 
are  required  (fig  3).  If  an  error  of  10~3  is  admissible 
then  only  13  control  points  are  needed.  These  errors  are 
measured  1.0  being  the  chord  of  the  airfoil.  The  results 
have  been  checked  by  using  different  solution  methods: 


1 3 5 7 9 11  13  15  17  19  21  23  25 


NB 


Figure  3:  Bezier  fitting  error  for  different  airfoils 


Gauss- Jordan,  Cholesky  and  QR. 

We  must  also  mention  that,  in  any  case,  given  to  the 
nature  of  matrix  [Ay],  there  are  no  additional  benefits 
in  using  more  than  30  control  points  due  to  numerical 
errors. 

As  a main  conclusion,  the  use  of  25  control  points  is  rec- 
ommended, which  represents  in  general  21  design  vari- 
ables by  taking  into  account  the  boundary  conditions  at 
the  leading  and  trailing  edges.  This  value  is  much  higher 
than  those  mostly  used  up  to  now,  but  we  consider  it  nec- 
essary for  the  solution  to  be  independent  of  the  design 
space. 


4 Flow  solver 

The  flow  solver  utilized  in  the  present  code  for  3-D 
computations  is  the  FL022vis  (INTA  version).  This 
code  employs  a viscous-inviscid  interaction  method  in 
a direct  iterative  way  to  compute  the  flow  past  swept 
wings  in  transonic  regime.  It  is  based  on  the  inviscid 
solver  FLO 2 2 by  Jameson  and  Caughey  [5]  and  the 
integral  boundary  layer  code  BL3D  by  Stock  [6],  recently 
modified  by  Yang  and  Wichmann  at  DLR  [7].  The 
inviscid  solver  FLO 2 2 has  been  modified  at  INTA  to 
render  the  trailing  edge  as  a grid  line  [8].  Additionally, 
the  coupling  procedure  and  the  smoothing  techniques 
for  the  displacement  thickness  have  been  improved.  A 
more  detailed  explanation  follows. 

4.1  Inviscid  solver 

The  inviscid  solver  is  a finite  differences  code  which  solves 
the  full  potential  equation  in  its  non- conservative  form 
in  a non  orthogonal  co-ordinate  system.  The  co-ordinate 
system  is  generated  by  a sequence  of  transformations.  In 
its  original  formulation,  the  spanwise  co-ordinate  lines 
were  aligned  with  the  leading  edge,  but  cut  across  the 
trailing  edge  in  a tapered  wing.  In  order  to  render  the 
trailing  edge  as  a grid  line,  an  additional  seeding  by  the 
local  wing  chord  has  been  introduced  in  each  of  the  wing 
sections.  In  the  resulting  co-ordinate  system,  both  the 
leading  and  trailing  edges  are  grid  lines.  New  stretch- 
ing transformations  have  been  used  outside  the  wing  to 
increase  the  distance  of  the  far  field  to  the  wing  and  to 
better  implement  the  boundary  conditions.  The  stretch- 
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Ala  DLR-F4,  CL=0.5,  Re=3*106 
Mach=0.75 


Figure  4:  RAE  2822  airfoil.  Cp  distribution  at 
M = 0.73,  Re  = 6.5  • 106  and  CL  = 0.8 


ing  in  normal  direction  at  transonic  flow  is  especially  im- 
portant. The  resultant  system  of  non  linear  algebraic 
equation  is  solved  by  an  iterative  relaxation  procedure 
SLOR  (Successive  Line  Over  Relaxation). 

4.2  Viscous  solver 

The  viscous  solver  is  an  integral  boundary  layer  code. 
It  solves  the  laminar  and  turbulent  boundary  layers  in  a 
general  non- orthogonal  curvilinear  coordinate  system.  In 
the  laminar  part,  a family  of  similarity  profiles  of  Falkner- 
Skan  is  used  in  the  streamwise  direction,  whilst  a com- 
bination of  these  profiles  are  used  in  crossflow  direction. 
For  the  turbulent  boundary  layer,  Coles  velocity  profiles 
of  two  parameters  are  used  in  the  streamwise  direction, 
and  Mager  or  Johnston  profiles  are  used  in  the  crossflow 
direction.  Transition  is  fixed  or  is  computed  by  a set 
of  empirical  correlations.  The  4th  — order  Runge-Kutta 
scheme  is  used  for  the  numerical  integration,  in  which 
the  ^-direction  is  the  marching  direction. 

4.3  Coupling  procedure 

The  concept  of  displacement  thickness  is  used  for  the  in- 
teraction of  the  external  flow  with  the  boundary  layer. 
This  method  only  provides  a proper  treatment  of  ‘weak 
interaction’  phenomena.  An  underrelaxation  technique  is 
introduced  to  avoid  instabilities,  and  a special  smooth- 
ing of  the  displacement  surface  by  means  of  cubic  Bezier 
splines  is  employed  to  avoid  strong  irregularities  in  the 
shock  region  and  in  the  trailing  edge. 

4.4  Validation 

The  flow  solver  has  been  widely  validated  after  its  im- 
provement. Two  sample  cases  in  transonic  flow  will  be 
shown  below.  It  is  important  to  compare  not  only  global 
values  but  also  pressure  distributions.  For  potential  flow, 
isentropic  flow  and  irrotational  flow  are  assumed.  Dis- 
crepancy with  the  correct  Rankine-Hugoniot  result  is 
small  if  the  shock  upstream  Mach  number  is  M\  < 1.3. 
Moreover,  the  viscous-inviscid  interaction  takes  into  ac- 
count only  weak  interaction  phenomena.  At  M\  >1.3 
incipient  separation  occurs,  thus  our  viscous-inviscid  cou- 
pling is  less  accurate. 


Figure  5:  DLR  F4  wing.  Cp  distributions  on  several 
span  wise  sections  at  M — 0.75,  Re  — 3 • 106  and 
CL  = 0.5 

4.4.1  RAE  2822  airfoil 

A limit  case  for  comparison  may  be  case  9 of  reference 
[9].  At  M = 0.73,  Re  = 6.5  • 106  and  a = 3.19  deg.  the 
boundary  layer  measurements  show  flow  close  to  separa- 
tion, and  the  shock  upstream  Mach  number  is  Mi  « 1.3. 
Figure  4 shows  the  comparison  of  experimental  and  com- 
puted pressure  distributions.  The  dashed  line  shows  the 
result  using  a grid  of  192-24*32  points,  and  the  solid  line 
represents  the  results  with  a mesh  of  192  -48-32  points. 
There  is  an  important  difference  in  the  shock  location. 
The  reason  is  that,  with  the  stretching  transformation 
used  in  the  normal  direction,  the  far  field  is  too  close  to 
the  wing  with  the  coarser  mesh.  A finer  grid  is  needed 
to  have  the  far  field  at  a sufficient  distance  to  avoid  any 
influence  of  the  far  field  into  the  supersonic  zone.  Al- 
though not  shown  here,  the  boundary  layer  parameters 
show  good  agreement  with  the  experimental  ones  before 
the  shock.  After  the  shock,  there  are  some  differences, 
especially  in  the  displacement  thickness. 

4.4.2  DLR-F4  wing 

Another  test  case  of  a typical  transonic  wing  will  be 
shown.  During  1997-98  INTA  has  participated  in  the 
A ERE  A F4  Model  Test  Programme  to  study  the  scale 
effects  on  the  DLR-F4  wing-body  model  in  the  new 
criogenic  european  tunnel  ETW  (European  Transonic 
Windtunnel).  The  range  of  Reynolds  number  was 
3 — 33  * 106  and  the  range  of  Mach  number  was  0.6  — 0.81. 
Thus,  the  amount  of  experimental  information  is  very 
important  and  permits  the  validation  and  assessment  of 
CFD  codes  [10,  11].  In  this  AEREA  campaign,  pressure 
measurements  were  not  made,  but  there  are  experimen- 
tal data  obtained  by  the  Garteur  Group  AG-01  [12]. 
Figine  5 shows  the  Cp  distribution  in  several  spanwise 
sections  of  the  DLR-F4  wing  at  M = 0.75,  Cl  = 0.5  and 
Re  = 3.0  • 106. 

The  computational  mesh  used  had  192  -48-64  grid 
points  and  the  computational  time  was  about  l|  hour 
in  a SGI  DN  10000  Power  Challenge  workstation  (300 
mfiops).  The  mesh  is  finer  than  in  the  previous  case 
to  take  into  account  the  large  variations  of  geometry 
in  spanwise  direction;  otherwise  the  shock  will  not  be 
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DLR-F4  wing,  C,  =0.5 
Influence  of  Reynolds  number 


Figure  6:  DLR-F4  wing.  Reynolds  effects  on  drag 
coefficient  at  Cj,  — 0.5 

properly  captured.  It  can  be  seen  that  the  comparison 
of  the  theoretical  results  with  the  experiments  is  good. 
There  are  some  differences  in  the  lower  side  due  to  the 
incapability  of  the  flow  solver  to  predict  the  influence  of 
the  body,  and  also  to  difficulties  in  the  implementation 
of  Kutta  condition. 

With  regard  to  global  values,  figure  6 shows  the  drag 
coefficient  at  M = 0.785  and  Cl  = 0.5  at  several 
Reynolds  numbers.  The  comparison  is  done  with  fixed 
transition  (5%/5%)  and  free  transition.  The  fuselage 
drag  has  been  computed  using  an  analytical  procedure. 
The  figure  shows  that  drag  is  overestimated  using  the 
code,  but  it  is  worth  noting  that  the  experimental 
results  have  not  been  corrected,  and  that  Reynolds 
number  effects  include  pseudo  Reynolds  number  effects 
due  to  aeroelastic  distorsion.  Two  curves  representing 
a translation  of  the  theoretical  curves  have  additionally 
been  plotted  in  the  figure  to  compare  the  trends  of  both 
theoretical  and  experimental  results.  It  can  be  observed 
that,  despite  the  results  at  Reynolds  number  of  3 • 106 
for  fixed  transition  and  3 — 6 • 106  for  free  transition,  the 
influence  of  Reynolds  number  in  drag  is  well  predicted. 
It  must  be  underlined  that  at  Reynolds  number  of  3 • 106 
a small  trailing  edge  separation  is  present  in  the  kink 
region  of  the  trailing  edge  [12],  thus  the  drag  increases. 

A general  conclusion  is  that  the  code  can  predict  the 
flow  past  wings  in  transonic  regime  at  attached  flow  con- 
ditions with  sufficient  accuracy.  These  are  the  cruising 
conditions.  The  advantage  of  using  this  kind  of  code  is 
that  an  optimization  code  with  this  flow  solver  will  not 
be  highly  time  consuming. 


5 Results 

To  show  the  capabilities  of  the  method  two  cases  (sub- 
sonic and  transonic)  have  been  tested.  But  before  any 
optimization  run  can  be  started  it  is  important  to  take 
into  account  several  considerations.  With  the  same  en- 
gineering objective  there  are  many  different  ways  to  pro- 
ceed. For  example,  the  aerodynamic  efficiency  could  be 
optimized  by  using  it  directly  as  objective  function,  or  by 
optimizing  Cl  with  Cd  as  constraint  or  viceversa.  But  it 
is  also  possible  that  not  all  of  those  three  strategies  will 
provide  a feasible  design,  if  any. 

On  the  other  hand,  when  viscous  effects  are  included  in 
the  analysis  the  changes  in  the  objective  function  can  be 


Figure  7:  Sensitivity  analysis  for  the  subsonic  case 


highly  non  linear  and  it  is  very  important  to  choose  the 
most  appropriate  values  for  the  initial  design  variables, 
not  so  small  that  the  objective  does  not  change,  and  not 
so  large  that  the  modified  geometry  can  not  be  analyzed. 
That’s  why  we  have  first  computed  the  sensitivities  of 
Cl  and  Cd  to  changes  in  the  design  variables.  It  also 
gives  us  an  idea  of  the  convergence  criteria  to  use  in  the 
optimization  run  in  order  to  avoid  as  many  analyses  as 
possible.  Lastly,  it  provides  useful  information  about 
the  convenience  or  not  for  scaling  the  variables  and  for 
applying  limiting  values. 

5.1  Subsonic  case 

The  first  exercise  has  been  the  optimization  of  the 
wing  of  an  RPV  called  SIVA  (Sistema  Integrado  de 
Vigilancia  Area)  designed  at  INTA.  The  flow  conditions 
are  M = 0.2  and  Re  = 1.5  • 106.  The  wing  has  no 
twist,  a surface  of  4.5  sq.m.,  5.81  m of  span,  0.905 
and  0.485  m of  chord  at  the  root  and  at  the  tip, 
respectively.  Among  a broad  set  of  different  airfoils,  the 
Eppler-580  was  selected.  At  cruise  conditions  the  an- 
gle of  attack  is  1.5  deg.  and  the  lift  coefficient  Cl  =0.56. 

By  using  the  nine  design  variables  described  in  § 3 (eqs. 
1 to  4)  the  sensitivities  computed  for  a change  of  10 “s 
in  each  variable  are  plotted  in  fig.  7.  There  is  an  slight 
decrease  in  lift  and  drag  when  the  thickness  increases, 
but  the  efficiency  in  general  increases,  and  more  so 

when  the  deformation  approaches  the  trailing  edge.  The 
opposite  happens  as  regards  the  camber  deformation. 

From  these  values  the  following  problem  has  been  posed: 
improvement  of  at  a weighted  rate  of  0.5  at  0.0 

deg.  and  0.5  at  2.0  deg.  with  the  constraint  of  having 
a section  area  not  smaller  than  that  of  the  original 
one.  By  running  a multipoint  optimization,  the  use  of 
an  arbitrary  constraint  in  the  leading  edge  radius  is 
avoided.  It  is  also  advisable  not  to  use  higher  values  of 
angle  of  attack  when  the  accuracy  of  the  results  gets 
worse.  The  flow  is  supposed  to  be  fully  turbulent. 

A design  with  the  airfoil  Eppler-580  as  starting  geometry 
was  carried  out.  A sequence  of  three  meshes,  from  coarse 
to  fine,  are  utilized  in  order  to  accelerate  the  convergence 
and  to  save  computing  time.  The  geometry  obtained 
is  plotted  in  fig.  8 compared  with  the  original  one.  As 
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Figure  8:  Geometries  of  both  original  and  designed 
wing  sections  for  the  subsonic  case 


Figure  9:  Efficiency  vs  lift  coefficient  for  the  origi- 
nal and  designed  wing  in  inviscid  and  viscous  flow 
(subsonic  case) 

expected  from  the  sensitivity  results  the  camber  has 
decreased  for  most  of  the  chord.  The  efficiencies  for  a 
range  of  Cl*s  are  plotted  in  fig.  9 where  the  results 
for  an  inviscid  run  are  also  included.  For  the  viscous 
computation  the  efficiency  is  improved  9.4%. 

The  pressure  distributions  for  the  mid  section  of  both 
wings  are  represented  in  fig.  10.  It  can  be  seen  that 
the  maximum  velocity  has  decreased,  which  is  usually 
positive  in  terms  of  drag.  Besides,  the  pressure  gradient 
around  the  leading  edge  on  the  lower  side  is  much 
smoother.  A further  study  should  be  made  about  the 
more  adverse  pressure  gradient  close  to  the  trailing  edge 
on  the  upper  side,  using  different  verification  methods. 

In  order  to  check  the  dependency  of  the  solution  on  the 
initial  geometry  an  additional  design  has  been  carried 
out  with  the  same  objective  and  constraints  as  before, 
but  starting  from  a conventional  wing  section  (NACA- 
0015),  The  results  are  plotted  in  figures  11  and  12.  It  can 
be  seen  that  even  in  this  case  the  aerodynamic  charac- 
teristics of  the  original  section  have  also  been  improved, 


Figure  10:  Pressure  coefficient  at  central  section  of 
both  original  and  designed  wing  (subsonic  case) 


Figure  11:  Geometries  of  NACA  0015,  original  and 
designed  wing  sections  for  the  subsonic  case 

showing  in  a broader  sense  the  capabilities  of  the  method. 

Lastly,  a test  case  using  Bezier  polynomia  has  been 
carried  out.  It  is  enough  using  13  control  points  (10 
design  variables)  for  an  error  of  10“3,  as  it  was  pointed 
out  in  figure  3.  A finer  grid  than  in  the  previous  cases 
has  been  used.  The  grid  is  a C-H  mesh  of  160  -24-32 
points,  fine  enough  in  a subsonic  case.  The  initial 
wing  section  is  again  the  NACA  0015  airfoil.  The 
total  number  of  design  cycles  needed  were  8,  with 
85  evaluations  of  the  objective  function.  The  search 
direction  was  that  of  the  BFGS  (quasi-Newton)  method. 
Additionally,  there  is  no  asumption  of  fully  turbulent 
flow.  The  transition  is  computed  by  means  of  empirical 
methods.  The  computing  time  was  24  hours  49  min 
in  a SGI  Power  Challenge  workstation  (300  mflops). 
Figure  13  shows  the  initial  section,  the  original  section 
of  the  SIVA  wing  and  the  designed  section  using  the 
Bezier  polynomia.  The  efficiencies  of  both  the  original 
and  designed  wing  are  plotted  in  figure  14.  In  order 
to  increase  confidence  in  the  new  design,  two  different 
codes  have  been  used  to  analyze  both  the  original  and 
designed  wings.  One  code  is  the  solver  employed  by  our 
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Figure  12:  Efficiency  vs  lift  coefficient  for  the  original 
and  designed  wing  with  NACA  0015  as  initial  section 
(subsonic  case) 


Figure  14:  Efficiency  vs  lift  coefficient  for  the  origi- 
nal and  designed  wing  using  Bezier  polynomia  with 
NACA  0015  as  initial  section  (subsonic  case) 


Figure  13:  Geometries  of  NACA  0015,  original  and 
designed  wing  sections  for  the  subsonic  case  using 
Bezier  polynomia 


optimization  method,  i.e.  FL022vis  code,  and  the  other 
one  is  a commercial  low  order  panel  method  coupled  to 
an  integral  boundary  layer  method  [VS AERO  code). 
This  code  computes  a lower  drag  than  the  former  one, 
but  it  is  worth  noting  that  both  codes  show  similar 
trends  for  both  wings.  There  are  some  differences  in 
the  high  Cl  values  range.  In  this  case  the  results  of 
the  finite  differences  code  are  less  accurate.  In  cruise 
condition  both  codes  predict  a similar  improvement  of 
the  efficiency. 

A final  comment  on  this  case  is  that,  a design  of  the  wing 
section  of  an  RPV  has  been  performed  starting  from  an 
initial  section  far  from  the  ‘optimal’  one.  The  result  is 
an  airfoil  with  similar  behaviour  of  the  selected  airfoil 
(Eppler-580)  at  cruise  conditions,  and  also  at  off- design 
conditions.  A very  few  number  of  constraints  have  been 
imposed.  This  result  shows  the  capabilities  of  the  opti- 
mization method  presented  in  this  paper.  The  cost  in 
terms  of  computational  effort  is  not  so  high. 


5.2  Transonic  case 

The  design  of  transonic  wings  is  more  time  consuming 
than  in  subsonic  flow.  In  order  to  simplify  the  opti- 
mization procedure,  the  design  shown  in  this  paper 
was  done  fixing  the  planform  and  twist,  and  the  wing 
sections  are  the  same  in  every  spanwise  section  of 
the  wing.  In  these  conditions,  the  improvement  of 
wing  features  will  depend  on  the  airfoil  features.  The 
wing  to  improve  is  a trapezoidal  wing  of  aspect  ratio 
AR  = 8.0,  leading  edge  sweep  angle  A = 27.5  deg., 
taper  ratio  A = 0.4  and  linear  twist.  A typical  grid 
for  wing  analysis  has  192  • 48  • 48  points.  For  the 
optimization,  a coarser  mesh  will  be  used,  because  of  the 
computational  time.  Figure  15  shows  the  transonic  wing. 

With  the  same  design  variables  as  used  in  the  subsonic 
case,  the  sensitivities  are  plotted  in  figs.  16, 17.  This  time 
comparing  the  inviscid  and  viscous  results. 

As  shown,  the  risk  of  obtaining  very  different  designs 
when  using  viscous  and  inviscid  analysis  is  high,  as 
the  sign  of  the  gradients  is  different  for  some  design 
variables  and  aerodynamic  coefficients.  For  instance, 
by  increasing  the  thickness  in  the  rear  part  of  the  wing 
section  the  drag  also  increases  when  the  analysis  is 
inviscid  but  decreases  for  a viscous  analysis. 

The  design  conditions  are  M = 0.82,  Re  = 6.0  * 106.  The 
goal  is  to  improve  (^^)  at  a weighted  rate  of  0.5  at  1.5 
deg.  and  0.5  at  3.0  deg.  with  the  constraint  of  having  a 
section  area  not  smaller  than  that  of  the  original  one. 

The  design  is  done  using  nine  design  variables  described 
in  § 3 (eqs.  1 to  4).  The  final  wing  section  can  be  seen 
in  fig.  18  along  with  the  original  one.  In  the  same  figure 
the  section  designed  by  running  the  analysis  as  inviscid 
is  plotted.  It  is  not  surprising  that  in  this  case  the  first 
impression  is  that  the  geometry  does  not  look  good 
due  to  the  solution  being  much  more  dependent  on  the 
design  variables  used  as  there  is  a greater  occurence 
of  local  minima.  The  need  to  take  viscous  effects  into 
account  is  clear,  especially  in  transonic  flow  where  the 
interaction  between  the  shock  wave  and  the  boundary 
layer  can  so  greatly  modify  the  designed  geometry. 

The  efficiencies  of  the  original  wing  and  the  viscous  de- 
sign are  plotted  in  fig.  19.  As  expected  in  any  tran- 
sonic case,  for  the  final  evaluation  of  both  two  wings  at 
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Figure  15:  Planform  view  of  the  transonic  wing. 
Typical  surface  mesh  for  analysis 


Figure  17:  Sensitivity  analysis  for  the  transonic  case 
(con’t) 


Figure  16:  Sensitivity  analysis  for  the  transonic  case 


Figure  18:  Geometries  of  both  original  and  designed 
wing  sections  for  the  transonic  case 


off-design  conditions  this  time  it  is  not  so  important  to 
include  viscous  effects  compared  with  the  subsonic  case 
(fig.  9)  as  the  relative  contribution  of  friction  drag  is 
smaller.  At  the  main  design  condition  the  efficiency  is 
improved  by  23.5%.  As  can  be  seen  in  fig. 20,  the  main 
changes  in  the  pressure  distribution  at  the  same  flow  con- 
dition are  located,  as  expected,  around  the  shock  wave, 
which  is  now  weaker,  while  it  remains  almost  the  same 
on  the  lower  surface. 


6 Conclusions 

A modular  optimization  code  for  wings  has  been 
described.  It  allows  any  combination  of  the  global 
aerodynamic  coefficients  to  be  used  as  objective  function 
along  with  a broad  set  of  geometrical  and  physical  con- 
straints. Even  though  it  is  possible  to  use  a prescribed 
pressure  distribution  as  an  objective,  it  is  not  necessary 
to  know  anything  about  it.  The  examples  shown  here 
have  been  run  on  a SGI  DN- 10000  Power  Challenge 
workstation  (300  mffops).  In  any  case,  most  designs  can 


be  obtained  from  one  day  to  another. 

There  are  many  more  capabilities  than  those  included 
in  this  paper,  especially  in  the  selection  of  the  number 
and  type  of  the  design  variables,  including  the  planform 
modification,  but  the  decision  will  be  provided  by  the 
application  to  more  different  cases.  At  the  moment,  the 
following  recommendations  can  be  mentioned: 

• From  our  experience  in  airfoil  design  (and  we  think 
it  can  also  be  extended  to  wings)  it  is  not  clear 
whether  in  any  of  the  cases  there  are  advantages  in 
using  high  level  codes  instead  of  low  level  ones,  tak- 
ing into  account  that  at  design  conditions  the  flow  is 
expected  to  be  attached  and  the  shock  waves,  if  any, 
are  weak.  Besides,  the  low  level  codes  are  generally 
more  robust. 

• On  the  contrary,  care  needs  to  be  taken  about  the 
size  of  the  design  space.  It  has  been  shown  that  un- 
less the  changes  in  the  geometry  are  expected  to  be 
small,  around  thirty  design  variables  for  a real  wing 
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Figure  19:  Efficiency  vs  lift  coefficient  for  the  orig- 
inal and  designed  wing  in  inviscid  and  viscous  flow 
(transonic  case) 


design  are  needed,  including  planform  modification. 

• There  is  no  engineering  design  without  the  possibil- 
ity of  running  multipoint  optimization. 

• Lastly,  even  though  the  viscous  contribution  can  be 
small,  in  any  case  the  analysis  should  include  vis- 
cous effects,  especially  in  the  transonic  cases  when 
no  pressure  distribution  is  prescribed,  even  if  the 
friction  drag  is  not  included  in  the  objective.  Oth- 
erwise, the  final  geometry  could  be  unfeasible  or  at 
least  very  dependent  of  the  design  variables  used. 

Finally,  in  the  short  term  we  plan  to  include  the  options 
already  available  in  the  2D  code  OPTPER  mentioned 
above  for  the  design  and  optimization  of  airfoils;  namely, 
to  add  the  option  of  free  transition  as  a new  design  vari- 
able [13]  and  to  use  control  theory  for  computing  the 
gradients  [14]. 
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DISCUSSION 


Session  III,  Paper  #21 


Dr  Nangia  (Nangia  Aero,  UK)  suggested  that  pitching  moment  should  be  included  in  the 
optimizations,  as  this  would  strongly  affect  Cl/Cd  conclusions  through  the  effect  on  trim  drag, 
etc... 


Mr  Monge  agreed,  pointing  out  that  the  option  to  take  account  of  pitching  moment  as  a 
constraint,  or  within  the  object  function,  is  already  included  in  his  code.  The  sample 
cases  presented  in  the  paper  were  only  intended  to  show  some  of  the  tool’s  capabilities. 


